rm(list = ls())
library(foreign)
library(sandwich)
library(CBPS)
library(MASS)
library(rms)
library(xtable)
library(plyr)
library(fmsb)
library(stargazer)
library(stats)
library(mice)
library(Amelia)
set.seed(786)

#Set working directory and load data
#setwd()
load('data.rData')

#Useful functions
checkBinaryTrait <- function(v, naVal="NA") {
	if (!is.numeric(v)) stop("Only numeric vectors are accepted.")
	vSet = unique(v)
	if (!missing(naVal)) vSet[vSet == naVal] = NA
	vSet = vSet[!is.na(vSet)]

	if (any(as.integer(vSet) != vSet)) FALSE
	else if (length(vSet) > 2) FALSE
	else TRUE
}
deci <- function(x, k, center = TRUE) {
	a <- format(round(x, k), nsmall=k)
	return(a)
}
decim <- function(x, k, center = TRUE) {
	a <- format(round(x, k), nsmall=k)
	if (center == TRUE & x >= 0.0){
		a <- as.character(a)
	}
	return(a)
}
testStatVector <- function(dat, balanceVars){
	returnVector <- NULL
	for (var in balanceVars){
		test.stat <- mean(dat[dat$cp == 1, var], na.rm = T) - mean(dat[dat$cp == 0, var], na.rm = T)
		returnVector <- c(returnVector, test.stat)
	}
	names(returnVector) <- balanceVars
	return(returnVector)
}
boot_lm <- function(dat, formul, cluster_var = 'ID', strata = NULL, R = 500){
	sample_coefs = NULL
	the_clusters <- as.character(unique(dat[[cluster_var]]))
	ref.model <- lm(formul, dat)
	coef_labels <- names(ref.model$coefficients)
	for (i in 1:R){
		boot_data <- data.frame()
		states.selected <- sample(x = the_clusters, size = length(the_clusters), replace = TRUE)
		bb <- table(states.selected)
		for(j in 1:max(bb)){
			selected.rows <- dat[[cluster_var]] %in% names(bb[bb %in% j])
			add_row <- dat[selected.rows, ]
			for(k in 1:j){
				boot_data <- rbind(boot_data, add_row)
			}
		}
		newrow <- rep(NA, length(coef_labels))
		names(newrow) <- coef_labels

		if(!is.null(strata)){
			boot_data$lm_weights <- NULL
			strat_levels <- levels(as.factor(boot_data[[strata]]))
			for (strat_level in strat_levels){
				this_strat <- boot_data[boot_data[[strata]] == strat_level,]
				pscore <- mean(this_strat$cp, na.rm = T)
				boot_data$lm_weights[boot_data[[strata]] == strat_level & boot_data$cp == 1] <- 1/pscore
				boot_data$lm_weights[boot_data[[strata]] == strat_level & boot_data$cp == 0] <- 1/(1-pscore)
			}
			#lm.out <- lm(formula = update(formul, . ~ . - factor(nAttempts)), data = boot_data, weights = lm_weights)
			lm.out <- lm(formula = formul, data = boot_data, weights = lm_weights)
			}else{lm.out <- lm(formula = formul, data = boot_data)}
		the.coefs <- summary(lm.out)$coefficients[, 1]
		for (thevar in names(the.coefs)){
			newrow[names(newrow) == thevar] <- the.coefs[names(the.coefs) == thevar]
		}
		sample_coefs <- rbind(sample_coefs, newrow)
	}
	rownames(sample_coefs) <- 1:R
	return(sample_coefs)
}
lm_mi <- function(dat_t, controls = NULL, reweigh = NULL, nAttempts_cntrl = TRUE, ctrend = F){
	coefs <- NULL
	ses <- NULL
	if (!is.null(controls)){
		reg.formul <- formula(paste('DV_diff ~ cp + factor(nAttempts)', ' + ', paste(controls, '.pt1', sep = '', collapse = ' + '), sep = ''))
	}else{
		reg.formul <- DV_diff ~ cp + factor(nAttempts)
	}
	if (nAttempts_cntrl == FALSE){
		reg.formul <-  update(reg.formul, . ~. -factor(nAttempts))
	}
	if (is.null(reweigh)){
		for (i in levels(dat_t$.imp)) {
		    in.imp <- dat_t$.imp == i
		    dt <- dat_t[in.imp, ]
		    did.out <- lm(reg.formul, data = dt, weights = NULL)
			new_coef <- summary(did.out)$coef[,1]
			new_se <- summary(did.out)$coef[,2]
			coefs <- rbind(coefs, new_coef)
			ses <- rbind(ses, new_se)
		}
	}else{
		for (i in levels(dat_t$.imp)) {
		    in.imp <- dat_t$.imp == i
		    dt <- dat_t[in.imp, ]
		    did.out <- lm(reg.formul, data = dt, weights = dt[[reweigh]])
			new_coef <- summary(did.out)$coef[,1]
			new_se <- summary(did.out)$coef[,2]
			coefs <- rbind(coefs, new_coef)
			ses <- rbind(ses, new_se)
		}
	}
	coefs_combined <- as.vector(mi.meld(coefs, ses)$q.mi)
	ses_combined <- as.vector(mi.meld(coefs, ses)$se.mi)
	zvals_t <- -as.vector(abs(coefs_combined/ses_combined))
	zTests <- 2*pnorm(zvals_t)

	a <- cbind(coefs_combined, ses_combined, 2*pnorm(zvals_t))
	rownames(a) <- colnames(mi.meld(coefs, ses)$q.mi)
	colnames(a) <- c('Coefs', 'Ses', 'PVal')
	return(list('coefs' = a, 'N' = nrow(dt), 'Nsucc' = sum(dt$cp), model = did.out))
}
vcovCluster <- function(model, cluster){
	if(nrow(model.matrix(model))!=length(cluster)){
	  stop("Cluster variable and model have different N")
	}
	M <- length(unique(cluster))
	N <- length(cluster)
	K <- model$rank
	dfc <- (M/(M - 1)) * ((N - 1)/(N - K ))
	uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
	rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
	return(rcse.cov)
}
panel_reg <- function(data.in, formul, ...){
	arguments <- list(...)
	if ('weights' %in% names(arguments)){
		if (!is.null(arguments[['weights']])){
			weights.in <- names(arguments)
			weights.lab <- unlist(arguments)
			data.in$weights <- data.in[[weights.lab]]
			did.out <- lm(formula = as.formula(formul), data = data.in, y = TRUE, weights = weights)
		}else{
			did.out <- lm(formula = as.formula(formul), data = data.in, y = TRUE, weights = NULL)
		}
	}else{
		did.out <- lm(formula = as.formula(formul), data = data.in, y = TRUE, weights = NULL)
	}
	#lm(as.formula(formul), data = data.in, y = TRUE, weights = data.in[[cus_weights]])
	robust_vcov <- vcovCluster(model = did.out, cluster = data.in$ccode)
	a <- coeftest(did.out, robust_vcov, df.residual(did.out))
	return(list(model = did.out, coefs = a))
}

#-------------------------------
#Variable description
#-------------------------------

#Variables ending with L1 are lagged 1 year
#Variables ending with L2 are lagged 2 years
#Variables ending with L3 are lagged 3 years, and so on
#cp denotes coup success
#cpFl denotes coup failure
#cpAt denotes coup attempt
#nAttempts counts the number of coup attempts
#e_regionpol denotes world region
#The .pt1 suffix denotes pre-treatment value of the variable

#-------------------------------
#TABLE 1
#-------------------------------

balanceVars <- c(
	'e_migdppclnL1', #GDP per capita
	'e_migdpgroL1', #GDP growth
	'ln.rexportspgdpL1', #Exports as share of GDP (log)
	'oil_binaryL1', #Oil and gas production
	'ethfracL1', #Ethnic fractionalization
	'ln.tpopL1', #Total population (logged)
	'tpopgroL1', #Population growth
	'cheibubDemoL1', #Democracy (CGV)
	'cheibubPresidL1', #Presidentialism (CGV)
	'cheibubMilL1',  #Military regime (CGV)
	'gwf_partyL1', #Party regime (GWF)
	'gwf_monarchL1', #Monarchy (GWF)
	'gwf_personalL1', #Personalist regime (GWF)
	'cheibubLegisL1', #Legislature (CGV)
	'ln.milperpcL1', #Per capita military personnel (log)
	'ln.milperL1', #Total military personnel (log)
	'ln.rCOWmilexL1', #Total military expenditures (log)
	'ln.CowmilexpgdpL1', #Military expenditures, share of GDP (log)
	'ln.rCOWmilexpcL1', #Per capita military expenditures (log)
	'ln.rCOWmilexpsoldL1', #Per soldier military expenditures (log)
	'cbd.changeCowMilexL1', #Change in military expenditures (cubic root)
	'effectivenumberL1', #Counterbalancing,  Pilster and Bohmelt
	'paramilitary.bscL1', #Counterbalancing, Belkin and Schofer
	'counterbal.pthL1', #Counterbalancing, Powell and Thyne
	'counterbalancingL1', #Counterbalancing, De Bruin
	'warL1', #Interstate war
	'CW', #Cold War
	'intra_warL1', #Intra-state war
	'prioL1', #UCDP Prio
	'OppNonViolBanksL1', #Unrest, non violent (Banks)
	'ncp.bsc', #Number of past successful coups
	'ncpFl.bsc', #Number of past failed coups
	'purgesL1', #Purges
	'latentmeanL1'# State repression
)

identify.binary.vars <- apply(dat[, balanceVars], 2, checkBinaryTrait)
bal_pvalues <- NULL
balanceTable <- matrix("", nrow = 0 * length(balanceVars), ncol = 5)
for (var in balanceVars){
	print (var)
	mean1 <- mean(dat[dat$cp == 1, var], na.rm = TRUE)
	mean0 <- mean(dat[dat$cp == 0, var], na.rm = TRUE)
	diff <- mean1 - mean0
	sd1 <- sd(dat[dat$cp == 1, var], na.rm = TRUE)
	sd0 <- sd(dat[dat$cp == 0, var], na.rm = TRUE)
	sddiff <- sqrt(sd1^2/nrow(dat[dat$cp == 0, ]) + sd0^2/nrow(dat[dat$cp == 0, ]))
	sd_pooled <- sd(dat[, var], na.rm = TRUE)
	if (identify.binary.vars[[var]] == FALSE){
		pvalue <- t.test(dat[dat$cp == 1, var], dat[dat$cp == 0, var]) $ p.value
	} else{
		pvalue <- chisq.test(table(dat[, c('cp', var)]), simulate.p.value = F) $ p.value
	}
	bal_pvalues <- c(bal_pvalues, pvalue)

	balanceTable <- rbind(balanceTable, c(
		var,
		decim(mean1, 3),
		decim(mean0, 3),
		decim(diff, 3),
		decim(pvalue, 3)
		)
	)
	balanceTable <- rbind(balanceTable, c(
		'',
		paste('(', decim(sd1, 3, center = FALSE), ')', sep = ''),
		paste('(', decim(sd0, 3, center = FALSE), ')', sep = ''),
		paste('(', decim(sddiff, 3, center = FALSE), ')', sep = ''),
		''#,
		#''
		)
	)
}

names(bal_pvalues) <- balanceVars
colnames(balanceTable) <- c('Variable', 'Success', 'Failure', 'Difference', 'p-val on difference')
Table1 <- print(
			 xtable(balanceTable,
			 caption = "Are Successful and Failed Coups Similar?",
			 align = 'llcccc',
			  label = 'balTab'),
		 include.rownames = FALSE,
		 caption.placement = 'top',
		 table.placement = "H",
		 #size = "\\scriptsize",
		 booktabs=TRUE,
		 floating = FALSE)

#Joint significance
prob_treatment <- sum(dat$cp)/nrow(dat)
monte_carlo_dat <- dat[, c(balanceVars)]

simulations <- 1:1000
Number_sign_at0.05 <- NULL
Number_sign_at0.01 <- NULL
N_greater_than.5 <- sum(bal_pvalues < .05)

for (simulation in simulations){
	bal_pvalues_sims <- NULL
	for (var in balanceVars){
		#print (var)
		treated <- sample(1:nrow(monte_carlo_dat), round(nrow(monte_carlo_dat)*prob_treatment))
		monte_carlo_dat$cp <- 0
		monte_carlo_dat$cp[treated] <- 1
		if( !(identify.binary.vars[names(identify.binary.vars) %in% var])){
			#t-test for continuous variables
			pvalue <- t.test(monte_carlo_dat[monte_carlo_dat$cp == 1, var], monte_carlo_dat[monte_carlo_dat$cp == 0, var]) $ p.value
		} else{
			#chi squared test for binary variables
			pvalue <- chisq.test(table(monte_carlo_dat[, c('cp', var)]), simulate.p.value = FALSE) $ p.value
		}
		bal_pvalues_sims <- c(bal_pvalues_sims, pvalue)
	}
	Number_sign_at0.05 <- c(Number_sign_at0.05, sum(bal_pvalues_sims < 0.05))
}
pLess.5 <- sum(Number_sign_at0.05 >= 3)/max(simulations)
pLess.5

#-------------------------------
#TABLE 2
#-------------------------------

#Generalized linear model with multiple imputation
glm_mi <- function(longDat, frml){
	analyses <- list()
	for (i in levels(longDat$.imp)) {
		if (as.integer(i) != 0){
			theDat <- longDat[longDat$.imp == i, ]
			analyses[[as.integer(i)]] <- with(theDat, glm(as.formula(frml), family = binomial))
			analyses[[as.integer(i)]] <- glm(as.formula(frml), family = binomial, data = theDat)
		}
	}
	return(as.mira(analyses))
}

balanceVars.glm <- c(
	'e_migdppclnL1',
	'e_migdpgroL1',
	'ln.rexportspgdpL1',
	'oil_binaryL1',
	'ethfracL1',
	'ln.tpopL1',
	'tpopgroL1',
	'cheibubDemoL1',
	'cheibubPresidL1',
	'cheibubMilL1',
	'cheibubLegisL1',
	'ln.milperpcL1',
	'ln.rCOWmilexpcL1',
	'counterbalancingL1',
	'warL1',
	'CW',
	'prioL1',
	'OppNonViolBanksL1',
	'ncp.bsc',
	'ncpFl.bsc',
	'purgesL1',
	'latentmeanL1'
	)

modelFit.null <- glm_mi(dat_mi, 'cp ~ 1')
modelFit.nullFE <- glm_mi(dat_mi, 'cp ~ e_regionpol')

form.short <- paste("cp ~ ", paste(c(balanceVars.glm[1:round(length(balanceVars.glm)/2)]), collapse = ' + '), collapse = "")
short.model <- glm_mi(longDat = dat_mi, frml = form.short)

form.shortFE <- paste("cp ~ ", paste(c('e_regionpol', balanceVars.glm[1:round(length(balanceVars.glm)/2)]), collapse = ' + '), collapse = "")
short.modelFE <-  glm_mi(dat_mi, form.shortFE)

form.long <- paste("cp ~ ", paste(c(balanceVars.glm), collapse = ' + '), collapse = "")
long.model <-  glm_mi(dat_mi, form.long)

form.longFE <- paste("cp ~ ", paste(c('e_regionpol', balanceVars.glm), collapse = ' + '), collapse = "")
long.modelFE <-  glm_mi(dat_mi, form.longFE)

cfs1 <- summary(pool(short.model))[, 'estimate']
cfs2 <- summary(pool(short.modelFE))[, 'estimate']
cfs3 <- summary(pool(long.model))[, 'estimate']
cfs4 <- summary(pool(long.modelFE))[, 'estimate']
names(cfs1) <- rownames(summary(pool(short.model)))
names(cfs2) <- rownames(summary(pool(short.modelFE)))
names(cfs3) <- rownames(summary(pool(long.model)))
names(cfs4) <- rownames(summary(pool(long.modelFE)))

ses1 <- summary(pool(short.model))[, 'std.error']
ses2 <- summary(pool(short.modelFE))[, 'std.error']
ses3 <- summary(pool(long.model))[, 'std.error']
ses4 <- summary(pool(long.modelFE))[, 'std.error']
names(ses1) <- rownames(summary(pool(short.model)))
names(ses2) <- rownames(summary(pool(short.modelFE)))
names(ses3) <- rownames(summary(pool(long.model)))
names(ses4) <- rownames(summary(pool(long.modelFE)))

pvls1 <- summary(pool(short.model))[, 'p.value']
pvls2 <- summary(pool(short.modelFE))[, 'p.value']
pvls3 <- summary(pool(long.model))[, 'p.value']
pvls4 <- summary(pool(long.modelFE))[, 'p.value']
names(pvls1) <- rownames(summary(pool(short.model)))
names(pvls2) <- rownames(summary(pool(short.modelFE)))
names(pvls3) <- rownames(summary(pool(long.model)))
names(pvls4) <- rownames(summary(pool(long.modelFE)))

cfs <- list(
	cfs1,
	cfs2,
	cfs3,
	cfs4
	)

se <- list(
	ses1,
	ses2,
	ses3,
	ses4
	)

pvals_tab2 <- c(
	deci(pool.compare(short.model, modelFit.null, method = 'likelihood', data = mice.out)$pvalue, 2),
	deci(pool.compare(short.modelFE, modelFit.null, method = 'likelihood', data = mice.out)$pvalue, 2),
	deci(pool.compare(long.model, modelFit.null, method = 'likelihood', data = mice.out)$pvalue, 2),
	deci(pool.compare(long.modelFE, modelFit.null, method = 'likelihood', data = mice.out)$pvalue, 2)
	)

pvalsFE_tab2 <- c(
	deci(pool.compare(short.model, modelFit.null, method = 'likelihood', data = mice.out)$pvalue, 2),
	deci(pool.compare(short.modelFE, modelFit.nullFE, method = 'likelihood', data = mice.out)$pvalue, 2),
	deci(pool.compare(long.model, modelFit.null, method = 'likelihood', data = mice.out)$pvalue, 2),
	deci(pool.compare(long.modelFE, modelFit.nullFE, method = 'likelihood', data = mice.out)$pvalue, 2)
	)

placeholderShort <- glm(as.formula(form.short), data = dat)
placeholderShortFE <- glm(as.formula(form.shortFE), data = dat)
placeholderLong <- glm(as.formula(form.long), data = dat)
placeholderLongFE <- glm(as.formula(form.longFE), data = dat)

Table2 <- stargazer(
			placeholderShort,
			placeholderShortFE,
			placeholderLong,
			placeholderLongFE,
			label = 'balReg',
			coef = cfs,
			se = se,
			dep.var.caption = "Dependent Variable: Coup Success",
			dep.var.labels = c('', '', '', ''),
			align = TRUE,
			omit = c('decade', 'region', 'Constant'),
			omit.stat = c('all'),
			title = 'Logistic Regression of Coup Success',
			add.lines = list(
				c(' Region Fixed Effects', 'No', 'Yes', 'No', 'Yes'),
				c(' Observations', deci(nrow(dat), 0), deci(nrow(dat), 0), deci(nrow(dat), 0), deci(nrow(dat), 0)),
				c(' $p$-value: $F$-test on all variables', pvals_tab2),
				c(' $p$-value: $F$-test on all variables (except FEs)', pvalsFE_tab2)
				),
			float = FALSE,
			notes.append = FALSE,
			notes = ''
		)

#--------------------------------
#Table 4
#--------------------------------

controls <- c(
				'cheibubMil',
				'ythblgap',
				'cheibubDemo',
				'ln.rCOWmilexpc',
				'CW',
				'counterbalancing',
				'ncp.bsc',
				'ncpFl.bsc'
			)
match.on <- c(
	 'e_migdppcln',
	 'cheibubLegis',
	 'ln.milperpc',
	 'counterbalancing',
	 'oil_binary',
	 'purges',
	 'ncp.bsc',
	 'ncpFl.bsc',
	 'OppNonViolBanks',
	 'ethfrac'
	 )

#Create weights
diff_dat$nAttempts[as.integer(as.character(diff_dat$nAttempts)) >= 3] <- 3
diff_temp1 <- subset(diff_dat, nAttempts ==1)
diff_temp2 <- subset(diff_dat, nAttempts ==2)
diff_temp3 <- subset(diff_dat, nAttempts ==3)
pscore1 <- mean(diff_temp1$cp)
pscore2 <- mean(diff_temp2$cp)
pscore3 <- mean(diff_temp3$cp)

diff_dat$lm_weights[as.integer(as.character(diff_dat$nAttempts)) == 1 & diff_dat$cp == 1] <- 1/pscore1
diff_dat$lm_weights[as.integer(as.character(diff_dat$nAttempts)) == 1 & diff_dat$cp == 0] <- 1/(1-pscore1)
diff_dat$lm_weights[as.integer(as.character(diff_dat$nAttempts)) == 2 & diff_dat$cp == 1] <- 1/pscore2
diff_dat$lm_weights[as.integer(as.character(diff_dat$nAttempts)) == 2 & diff_dat$cp == 0] <- 1/(1-pscore2)
diff_dat$lm_weights[as.integer(as.character(diff_dat$nAttempts)) == 3 & diff_dat$cp == 1] <- 1/pscore3
diff_dat$lm_weights[as.integer(as.character(diff_dat$nAttempts)) == 3 & diff_dat$cp == 0] <- 1/(1-pscore3)

diff_dat$ID <- 1:nrow(diff_dat)

#Covariance balancing propensity score
diff_dat_mi$ID <- 1:nrow(diff_dat_mi)
diff_dat_mi$ccode <- factor(as.character(diff_dat_mi$ccode))
diff_dat_mi$w <- NA
match.formul <- formula(paste('cp ~ ', paste(paste(match.on, '.pt1', sep = '', collapse = ' + '), sep = ''), ' + nAttempts', sep = ''))
for (i in levels(diff_dat_mi$.imp)) {
	in.imp <- diff_dat_mi$.imp == i
	m.out <- CBPS(match.formul, standardize = F, ATT = 0, data = diff_dat_mi[in.imp,])
	diff_dat_mi$w[in.imp] <- m.out$weights
	diff_dat_mi$ID[in.imp] <- diff_dat$ID
	diff_dat_mi$lm_weights[in.imp] <- diff_dat$lm_weights

}

#Baseline
#------------------------------------------
library(lmtest)
formul <- DV_diff ~ cp + factor(nAttempts)
did.out_boot <- boot_lm(dat = diff_dat, formul = formul, cluster_var = 'ID', strata = 'nAttempts')
did.out <- lm(formul, data = diff_dat, weights = lm_weights)
coefs.tab4.out <- coef(did.out)
SEs.tab4.out <- apply(did.out_boot, 2, sd)
pvalues.tab4.out <- 2*pnorm(-abs(coefs.tab4.out/SEs.tab4.out))

estims.cov <- lm_mi(dat_t = diff_dat_mi, controls = controls, reweigh = 'lm_weights')
coefs.tab4.cov.out <- estims.cov$coefs[,1]
SEs.tab4.cov.out <- estims.cov$coefs[,2]
pvalues.tab4.cov.out <- estims.cov$coefs[,3]
N.cov.out <- estims.cov$N
Nsucc.cov.out <- estims.cov$Nsucc

#Matching
#------------------------------------------
estims.match <- lm_mi(dat_t = diff_dat_mi, controls = NULL, reweigh = 'w')
coefs.tab4.match.out <- estims.match$coefs[,1]
SEs.tab4.match.out <- estims.match$coefs[,2]
pvalues.tab4.match.out <- estims.match$coefs[,3]
N.match.out <- estims.match$N
Nsucc.match.out <- estims.match$Nsucc

estims.match.cov <- lm_mi(dat_t = diff_dat_mi, controls = controls, reweigh = 'w')
coefs.tab4.match.cov.out <- estims.match.cov$coefs[,1]
SEs.tab4.match.cov.out <- estims.match.cov$coefs[,2]
pvalues.tab4.match.cov.out <- estims.match.cov$coefs[,3]
N.match.cov.out <- estims.match.cov$N
Nsucc.match.cov.out <- estims.match.cov$Nsucc

#Trends
#------------------------------------------
formul_trend <- DV_diff ~ cp + factor(nAttempts) + trend.pt1
did.out_boot_trend <- boot_lm(dat = diff_dat, formul = formul_trend, cluster_var = 'ID', strata = 'nAttempts')
did.out_trend <- lm(formul_trend, data = diff_dat, weights = lm_weights)
coefs.tab4.trends.out <- coef(did.out_trend)
SEs.tab4.trends.out <- apply(did.out_boot_trend, 2, sd)
pvalues.tab4.trends.out <- 2*pnorm(-abs(coefs.tab4.trends.out/SEs.tab4.trends.out))
N.trend.out <- nrow(diff_dat)
Nsucc.cov.trend.out <- sum(diff_dat$cp)

estims.trendsqrd <- lm_mi(dat_t = diff_dat_mi, controls = c(controls, 'trend', 'trendsqrd'), reweigh = 'lm_weights')
coefs.tab4.trendsSqrd.out <- estims.trendsqrd$coefs[,1]
SEs.tab4.trendsSqrd.out <- estims.trendsqrd$coefs[,2]
pvalues.tab4.trendsSqrd.out <- estims.trendsqrd$coefs[,3]
N.trendssqr.out <- estims.trendsqrd$N
Nsucc.trendssqr.out <- estims.trendsqrd$Nsucc

coefsList <- list(
	coefs.tab4.out,
	coefs.tab4.cov.out,
	coefs.tab4.match.out,
	coefs.tab4.match.cov.out,
	coefs.tab4.trends.out,
	coefs.tab4.trendsSqrd.out
	)

seList <- list(
	SEs.tab4.out,
	SEs.tab4.cov.out,
	SEs.tab4.match.out,
	SEs.tab4.match.cov.out,
	SEs.tab4.trends.out,
	SEs.tab4.trendsSqrd.out
	)

pList <- list(
	pvalues.tab4.out,
	pvalues.tab4.cov.out,
	pvalues.tab4.match.out,
	pvalues.tab4.match.cov.out,
	pvalues.tab4.trends.out,
	pvalues.tab4.trendsSqrd.out
	)

successes <- c(
	sum(diff_dat$cp),
	Nsucc.cov.out,
	Nsucc.match.out,
	Nsucc.match.cov.out,
	Nsucc.cov.trend.out,
	Nsucc.trendssqr.out
	)

attempts <- c(
	nrow(diff_dat),
	N.cov.out,
	N.match.out,
	N.match.cov.out,
	N.trend.out,
	N.trendssqr.out
	)

DiD_state <- stargazer(
	did.out,
	estims.cov$model,
	estims.match$model,
	estims.match.cov$model,
	did.out_trend,
	estims.trendsqrd$model,
	coef = coefsList,
	se = seList,
	p = pList,
	p.auto = FALSE,
	omit = c('factor', 'region', 'Constant'),
	dep.var.caption  = "\\emph{Dependent variable}: First Difference in State Repression",
	dep.var.label  = c('', '', '', '', '', ''),
	column.labels   = c('Baseline Models', "With CBPS Weighing", 'With Trends'),
	column.separate = c(2, 2, 2),
	dep.var.labels.include  = FALSE,
	align = TRUE,
	omit.stat = c('rsq', 'f', 'ser', 'adj.rsq', 'n'),
	add.lines = list(
		c("CBPS Weighting", '$No$', '$No$', '$Yes$', '$Yes$', '$No$', '$No$'),#, '$No$', '$No$', '$Yes$', '$Yes$'),
		c("Observations", attempts),#, '$No$', '$No$', '$Yes$', '$Yes$'),
		c("Successful Coups", successes)#, '$No$', '$No$', '$Yes$', '$Yes$'),
		),
	notes.append = FALSE,
	float = FALSE
	)

#------------------------------
#Figure 1
#------------------------------

dat_AR <- subset(diff_dat, cheibubDemo.pt1 == 0)
dat_AR$ID <- 1:nrow(dat_AR)

dat_AR_n1 <- subset(dat_AR, nAttempts == 1)
dat_AR_n2 <- subset(dat_AR, nAttempts == 2)
dat_AR_n3 <- subset(dat_AR, nAttempts == 3)
pscore1_AR <- mean(dat_AR_n1$cp)
pscore2_AR <- mean(dat_AR_n2$cp)
pscore3_AR <- mean(dat_AR_n3$cp)

dat_AR$lm_weights[dat_AR$nAttempts == 1 & dat_AR$cp == 1] <- 1/pscore1_AR
dat_AR$lm_weights[dat_AR$nAttempts == 1 & dat_AR$cp == 0] <- 1/(1-pscore1_AR)
dat_AR$lm_weights[dat_AR$nAttempts == 2 & dat_AR$cp == 1] <- 1/pscore2_AR
dat_AR$lm_weights[dat_AR$nAttempts == 2 & dat_AR$cp == 0] <- 1/(1-pscore2_AR)
dat_AR$lm_weights[dat_AR$nAttempts == 3 & dat_AR$cp == 1] <- 1/pscore3_AR
dat_AR$lm_weights[dat_AR$nAttempts == 3 & dat_AR$cp == 0] <- 1/(1-pscore3_AR)

did.out_boot_AR <- boot_lm(dat = dat_AR, formul = DV_diff ~ cp + factor(nAttempts), cluster_var = 'ID', strata = 'nAttempts')
l90_autoc <- quantile(did.out_boot_AR[,colnames(did.out_boot_AR)=='cp'], 0.05)
u90_autoc <- quantile(did.out_boot_AR[,colnames(did.out_boot_AR)=='cp'], 0.95)

did.out.AR <- lm(DV_diff ~ cp + factor(nAttempts),  dat_AR, weights = lm_weights)
coef_AR <- coef(did.out.AR)[names(coef(did.out.AR)) == 'cp']

dat_DEM <- subset(diff_dat, cheibubDemo.pt1 == 1)
dat_DEM$nAttempts[dat_DEM$nAttempts == 3] <- 2

dat_DEM_n1 <- subset(dat_DEM, nAttempts == 1)
dat_DEM_n2 <- subset(dat_DEM, nAttempts == 2)
pscore1_DEM <- mean(dat_DEM_n1$cp)
pscore2_DEM <- mean(dat_DEM_n2$cp)

dat_DEM$lm_weights[dat_DEM$nAttempts == 1 & dat_DEM$cp == 1] <- 1/pscore1_DEM
dat_DEM$lm_weights[dat_DEM$nAttempts == 1 & dat_DEM$cp == 0] <- 1/(1-pscore1_DEM)
dat_DEM$lm_weights[dat_DEM$nAttempts == 2 & dat_DEM$cp == 1] <- 1/pscore2_DEM
dat_DEM$lm_weights[dat_DEM$nAttempts == 2 & dat_DEM$cp == 0] <- 1/(1-pscore2_DEM)

did.out_boot_DEM <- boot_lm(dat = dat_DEM, formul = DV_diff ~ cp + factor(nAttempts), cluster_var = 'ID', strata = 'nAttempts')
l90_demo <- quantile(did.out_boot_DEM[,colnames(did.out_boot_DEM)=='cp'], 0.05)
u90_demo <- quantile(did.out_boot_DEM[,colnames(did.out_boot_DEM)=='cp'], 0.95)

did.out.DEM <- lm(DV_diff ~ cp + factor(nAttempts),  dat_DEM, weights = lm_weights)
coef_DEM <- coef(did.out.DEM)[names(coef(did.out.DEM)) == 'cp']

coefs <- c(coef_AR, coef_DEM)
l90 <- c(l90_autoc, l90_demo)
u90 <- c(u90_autoc, u90_demo)

pdf('Figure1.pdf', width = 6, height = 4)
	par(mgp=c(2.2,0.45,0), tcl=-0.4, mar=c(1.9,3.6,1.1,1.1))
	plot(x = c(1, 2), y = coefs, xlim = c(0.5, 2.5), ylim = c(-0.1, 0.35), xaxt = "n", yaxt = "n", ylab = 'Effect of Coup Success', las = 2, xlab = '')
	axis(side = 1, at = c(1, 2), labels = c('Coups in Autocracies', 'Coups in Democracies'))
	axis(side = 2, at = c(0.4, 0.3, 0.2, 0.1, 0, -0.1, -0.2), labels = c(0.4, 0.3, 0.2, 0.1, 0, -0.1, -0.2), las = 3)
	abline(h = 0, lty = 2, lwd = 0.5)
	arrows(x0 = c(1, 2), y0 = l90, y1 = u90, angle = 90, code = 3)
dev.off()

#------------------------------
#Figure 2
#------------------------------
library(mice)

analyses <- list()
for (imp in seq(1, 10)){
	diff_dat_t <- diff_dat
	diff_dat_t$latent.pt1 <- diff_dat_t[[paste('posterior_latentmean', imp, sep = '')]]
	did.out.marg.ef.dat <- na.omit(diff_dat_t[, c('DV_diff', 'cp', 'nAttempts', 'ccode', 'latent.pt1', 'Year', 'lm_weights')])
	did.out.marg.ef <- lm(DV_diff ~ cp*(latent.pt1) + factor(nAttempts), data = did.out.marg.ef.dat, weights = lm_weights)
	analyses[[imp]] <- did.out.marg.ef
}

pld <- pool(analyses)$pooled
rnms1 <- which(rownames(pld) == 'cp')
rnms2 <- which(rownames(pld) == 'cp:latent.pt1')
pld_cp_est <- pld[c(rnms1, rnms2), 'estimate']
pld_cp_var <- pld[c(rnms1, rnms2), 't']
pool_varVcov <- matrix(nrow = 2, ncol = 2, 0)
diag(pool_varVcov) <- (pld_cp_var)
draws <- mvrnorm(n = 100000, mu = pld_cp_est, Sigma = pool_varVcov)
colnames(draws) <- c('cp', 'cp:latent.pt1')

latentmean_levels <-  seq(-4, 4, 0.01)
frstTerm <- as.matrix(data.frame(rep(1, length(latentmean_levels)), latentmean_levels))
marginal_effects <- frstTerm %*% t(draws)
means <- apply(marginal_effects, 1, mean)
l90 <- apply(marginal_effects, 1, function(x) quantile(x, 0.05))
u90 <- apply(marginal_effects, 1, function(x) quantile(x, 0.95))
l95 <- apply(marginal_effects, 1, function(x) quantile(x, 0.025))
u95 <- apply(marginal_effects, 1, function(x) quantile(x, 0.975))

pdf('Figure2.pdf', width = 16*0.5, height = 12*0.3)
	par(mgp=c(2, 0.5, 0), tcl=-0.4, mar=c(3.1,3.6,1.1,1.1))
	attach(mtcars)
	plot(x = latentmean_levels,
		y = means,
		type = 'l',
		xlim = c(-2, 3),
		ylim = c(-0.15, 0.33),
		xlab = 'Pre-Treatment Level of State Repression',
		ylab = 'Effect of Coup Success',
		xaxt = 'n',
		yaxt = 'n',
		lty = 1,
		lwd = 1
		)
	axis(1, at = -2:3, labels = -2:3)
	axis(2, at = seq(-0.2, 0.3, by = 0.1), labels = seq(-0.2, 0.3, by = 0.1), las = 3)
	polygon(c(latentmean_levels, rev(latentmean_levels)),
		c((l90), rev((u90))),
		border = FALSE,
		fillOddEven = TRUE,
		col = rgb(171/255, 171/255, 171/255, alpha= 0.55))
	polygon(c(latentmean_levels, rev(latentmean_levels)),
		c((l95), rev((u95))),
		border = FALSE,
		fillOddEven = TRUE,
		col = rgb(171/255, 171/255, 171/255, alpha= 0.25))
	abline(0, 0, lty =1, lwd = 0.5)
	rug(diff_dat$latentmean.pt1)
	abline(v = 0, lty = 2, col = 'grey')
	abline(v = 1, lty = 2, col = 'grey')
	abline(v = 2, lty = 2, col = 'grey')
	abline(v = 3, lty = 2, col = 'grey')
	abline(v = -1, lty = 2, col = 'grey')
	abline(v = -2, lty = 2, col = 'grey')
	abline(v = -3, lty = 2, col = 'grey')
dev.off()


#------------------------------
#Figure 3
#------------------------------
t.min <- -3
t.max <- 2

#nYears counts the number of year that an observation lasts (in the differenced datasets)
pre1_suc <- mean(diff_dat$latentL1[diff_dat$cp == 1 & diff_dat$nYears == 1], na.rm = T)
pre1_fal <- mean(diff_dat$latentL1[diff_dat$cp == 0 & diff_dat$nYears == 1], na.rm = T)

pre2_suc <- mean(diff_dat$latentL2[diff_dat$cp == 1 & diff_dat$nYears == 1], na.rm = T)
pre2_fal <- mean(diff_dat$latentL2[diff_dat$cp == 0 & diff_dat$nYears == 1], na.rm = T)

pre3_suc <- mean(diff_dat$latentL3[diff_dat$cp == 1 & diff_dat$nYears == 1], na.rm = T)
pre3_fal <- mean(diff_dat$latentL3[diff_dat$cp == 0 & diff_dat$nYears == 1], na.rm = T)

post1_suc <- mean(diff_dat$latentLead1[diff_dat$cp == 1 & diff_dat$nYears == 1], na.rm = T)
post1_fal <- mean(diff_dat$latentLead1[diff_dat$cp == 0 & diff_dat$nYears == 1], na.rm = T)

post2_suc <- mean(diff_dat$latentLead2[diff_dat$cp == 1 & diff_dat$nYears == 1], na.rm = T)
post2_fal <- mean(diff_dat$latentLead2[diff_dat$cp == 0 & diff_dat$nYears == 1], na.rm = T)

suc0 <- mean(diff_dat$dv[diff_dat$cp == 1 & diff_dat$nYears == 1], na.rm = T)
fal0 <- mean(diff_dat$dv[diff_dat$cp == 0 & diff_dat$nYears == 1], na.rm = T)

meanSuccL <- c(pre3_suc, pre2_suc, pre1_suc)
meanSuccLead <- c(post1_suc, post2_suc)

meanFailL <- c(pre3_fal, pre2_fal, pre1_fal)
meanFailLead <- c(post1_fal, post2_fal)

successes <- c(meanSuccL, mean0 = suc0, meanSuccLead)
failures <- c(meanFailL, mean0 = fal0, meanFailLead)

dt <- data.frame(suc = successes, fail = failures, time = t.min:t.max)

newdata1 <- data.frame(
	time = t.min:0
	)
newdata2 <- data.frame(
	time = 0:t.max
	)

succes_pre_lm <- predict(lm(suc ~ time, data = dt, subset =  time<0), newdata = newdata1)
succes_pst_lm <- predict(lm(suc ~ time, data = dt, subset =  time>0), newdata = newdata2)
failure_pre_lm <-predict(lm(fail ~ time, data = dt, subset =  time<0), newdata = newdata1)
failure_pst_lm <-predict(lm(fail ~ time, data = dt, subset =  time>0), newdata = newdata2)

pdf('Figure3.pdf', width = 9, height = 4)
	plot(x = t.min:t.max, y = successes, xlim = c(t.min, t.max), pch = 15 , ylim = c(0.4, 0.79), xaxt = "n", yaxt = "n", xlab = '', ylab = '')
	lines(x = t.min:t.max, y = failures, type = 'p', pch = 20, col = 'grey')

	lines(x = t.min:0, y = succes_pre_lm)
	lines(x = t.min:0, y = failure_pre_lm, col = 'grey')
	lines(x = 0:t.max, y = succes_pst_lm)
	lines(x = 0:t.max, y = failure_pst_lm, col = 'grey')

	abline(v = 0, lty=2, col = 'grey')
	legend(x= -2, y = 1.2, c('Successful Coups', 'Failed Coups'), col = c('black', 'grey'), pch = c(15, 20), lty = c(1, 4), cex=  .66)
	mtext(side = 1, text = "Years Since or Until Coup Attempt", line = 1.9, cex = .86)
	mtext(side = 2, text = "Average State Repression", line = 2.5, cex = .86, las = 3)
	axis(1, at = seq(t.min, t.max, by = 1), labels = seq(t.min, t.max, by = 1), las = 1)
	axis(2, at = seq(-2, 2, by = 0.2), labels = seq(-2, 2, by = 0.2), las = 2)
dev.off()

#------------------------------
#Table 5
#------------------------------

D_bsc <- read.csv('ts_bsc.csv')
bsc.keep <- !is.na(D_bsc$cpAtL1)&!is.na(D_bsc$cpAt)&!is.na(D_bsc$cpAtLead1)&!is.na(D_bsc$latentmean)
formul.A <- latentmean ~ cpAtLead1 + factor(Year) + trend*factor(ccode) + trendsqrd*factor(ccode) + trendcube*factor(ccode) # + trendfrth*factor(ccode) + shrYth1524
formul.B <- latentmean ~ cpAtL1 + factor(Year) + trend*factor(ccode) + trendsqrd*factor(ccode) + trendcube*factor(ccode) # + trendfrth*factor(ccode) + shrYth1524

ts.bsc_t <- D_bsc[bsc.keep, ]
ts.bsc_t$ccode <- factor(ts.bsc_t$ccode)
lm.A.bsc <- panel_reg(dat = ts.bsc_t, formul = formul.A, weights = NULL)
lm.B.bsc <- panel_reg(dat = ts.bsc_t, formul = formul.B, weights = NULL)
cfs.A.bsc <- lm.A.bsc$coefs[, 1:4]
cfs.B.bsc <- lm.B.bsc$coefs[, 1:4]

cfs.A.bsc[grepl('^cp', rownames(cfs.A.bsc)), ]
cfs.B.bsc[grepl('^cp', rownames(cfs.B.bsc)), ]

coef.fe <- list(
	cfs.A.bsc[,1],
	cfs.B.bsc[,1]
	)

se.fe <- list(
	cfs.A.bsc[,2],
	cfs.B.bsc[,2]
	)

pval.fe <- list(
	cfs.A.bsc[,4],
	cfs.B.bsc[,4]
	)

Table5 <- stargazer(
		lm.A.bsc$model,
		lm.B.bsc$model,
		se = se.fe,
		p = pval.fe,
		omit = c('ccode', 'Year', 'Constant', 'trend'),
		dep.var.caption  = "Dependent Variable: State Repression",
		column.separate = c(2, 2),
		dep.var.labels  = c(''),
		align = TRUE,
		title = 'Fixed-Effects Regression of Coup Attempt with Leads and Lags',
		order = c( 'cp.Lead3', 'cp.Lead2', 'cp.Lead1', 'cpAt', 'cpAt.L1', 'cpAt.L2', 'cpAt.L3'),
		omit.stat = c('rsq', 'f', 'ser'),
		add.lines = list(
			c('Country FEs', '$Yes$', '$Yes$'),#, '$Yes$', '$Yes$'),
			c('Year FEs', '$Yes$', '$Yes$'),#, '$Yes$', '$Yes$'),
			c('Cubic Trends by Country', '$Yes$', '$Yes$')#,#, '$Yes$', '$Yes$'),
			),
		notes.append = FALSE,
		label = 'did.placeb.bsc',
		float = FALSE
		)


#------------------------------
#Table 6
#------------------------------

D_pth <- read.csv('ts_pth.csv')
pth.keep <- !is.na(D_pth$cpAtL1)&!is.na(D_pth$cpAt)&!is.na(D_pth$cpAtLead1)&!is.na(D_pth$latentmean)
dat_pth <- D_pth[pth.keep, ]
dat_pth$ccode <- factor(dat_pth$ccode)

dat_pth.CW1 <- subset(dat_pth, Year <= 1991)
dat_pth.CW1$ccode <- factor(dat_pth.CW1$ccode)
dat_pth.CW0 <- subset(dat_pth, Year > 1991)
dat_pth.CW0$ccode <- factor(dat_pth.CW0$ccode)

lm.A.pth.CW1 <- panel_reg(dat_pth.CW1, formul.A)
lm.B.pth.CW1 <- panel_reg(dat_pth.CW1, formul.B)
lm.A.pth.CW0 <- panel_reg(dat_pth.CW0, formul.A)
lm.B.pth.CW0 <- panel_reg(dat_pth.CW0, formul.B)

cfs.A.pth.CW1 <- lm.A.pth.CW1$coefs[, 1:4]
cfs.B.pth.CW1 <- lm.B.pth.CW1$coefs[, 1:4]
cfs.A.pth.CW0 <- lm.A.pth.CW0$coefs[, 1:4]
cfs.B.pth.CW0 <- lm.B.pth.CW0$coefs[, 1:4]

cfs.A.pth.CW1[grepl('^cp', rownames(cfs.A.pth.CW1)), ]
cfs.B.pth.CW1[grepl('^cp', rownames(cfs.B.pth.CW1)), ]
cfs.A.pth.CW0[grepl('^cp', rownames(cfs.A.pth.CW0)), ]
cfs.B.pth.CW0[grepl('^cp', rownames(cfs.B.pth.CW0)), ]

coef.fe <- list(
	cfs.A.pth.CW1[,1],
	cfs.B.pth.CW1[,1],
	cfs.A.pth.CW0[,1],
	cfs.B.pth.CW0[,1]
	)

se.fe <- list(
	cfs.A.pth.CW1[,2],
	cfs.B.pth.CW1[,2],
	cfs.A.pth.CW0[,2],
	cfs.B.pth.CW0[,2]
	)

pval.fe <- list(
	cfs.A.pth.CW1[,4],
	cfs.B.pth.CW1[,4],
	cfs.A.pth.CW0[,4],
	cfs.B.pth.CW0[,4]
	)

DiD_placeb.pth_t <- stargazer(
		lm.A.pth.CW1$model,
		lm.B.pth.CW1$model,
		lm.A.pth.CW0$model,
		lm.B.pth.CW0$model,
		se = se.fe,
		p = pval.fe,
		omit = c('ccode', 'Year', 'Constant', 'trend'),
		dep.var.caption  = "Dependent Variable: State Repression",
		column.labels   = c('During the Cold War', 'After the Cold War'),
		column.separate = c(2, 2),
		dep.var.labels  = c(''),
		align = TRUE,
		title = 'Fixed-Effects Regression of Coup Attempt during and after the Cold War, Powell and Thyne Data',
		order = c( 'cp.Lead3', 'cp.Lead2', 'cp.Lead1', 'cpAt', 'cpAt.L1', 'cpAt.L2', 'cpAt.L3'),
		omit.stat = c('rsq', 'f', 'ser'),
		add.lines = list(
			c('Country FEs', '$Yes$', '$Yes$', '$Yes$', '$Yes$'),
			c('Year FEs', '$Yes$', '$Yes$', '$Yes$', '$Yes$'),
			c('Cubic Trends by Country', '$Yes$', '$Yes$', '$Yes$', '$Yes$')
			),
		notes.append = FALSE,
		float = FALSE
		)

#------------------------------
#Table 6
#------------------------------

formul <- latentmean ~ cpLead2 + cpLead1 + cp + cpL1 + cpL2 + cpL3 + cpL4 + cpL5 + factor(ccode) + trend*factor(ccode) +factor(Year)+ trendsqrd*factor(ccode) + trendcube*factor(ccode)
minYear <- 1950; maxYear <- 2000
placeb.dat.bsc <- na.omit(D_bsc[D_bsc$Year >= minYear & D_bsc$Year <= maxYear, c('ccode', all.vars(formul))])
placeb.dat.bsc$ccode <- factor(placeb.dat.bsc$ccode)
placeb.reg.bsc <- lm(formul, data = placeb.dat.bsc)
CovVar.bsc <- vcovCluster(model = placeb.reg.bsc, cluster = placeb.dat.bsc$ccode)
results.bsc <- coeftest(placeb.reg.bsc, CovVar.bsc, df.residual(placeb.reg.bsc))

Means.bsc <- c(
	results.bsc[,1][names(results.bsc[,1]) == 'cpLead2'],
	results.bsc[,1][names(results.bsc[,1]) == 'cpLead1'],
	results.bsc[,1][names(results.bsc[,1]) == 'cp'],
	results.bsc[,1][names(results.bsc[,1]) == 'cpL1'],
	results.bsc[,1][names(results.bsc[,1]) == 'cpL2'],
	results.bsc[,1][names(results.bsc[,1]) == 'cpL3'],
	results.bsc[,1][names(results.bsc[,1]) == 'cpL4'],
	results.bsc[,1][names(results.bsc[,1]) == 'cpL5']
	)

SEs.bsc <- c(
	results.bsc[,2][names(results.bsc[,2]) == 'cpLead2'],
	results.bsc[,2][names(results.bsc[,2]) == 'cpLead1'],
	results.bsc[,2][names(results.bsc[,2]) == 'cp'],
	results.bsc[,2][names(results.bsc[,2]) == 'cpL1'],
	results.bsc[,2][names(results.bsc[,2]) == 'cpL2'],
	results.bsc[,2][names(results.bsc[,2]) == 'cpL3'],
	results.bsc[,2][names(results.bsc[,2]) == 'cpL4'],
	results.bsc[,2][names(results.bsc[,2]) == 'cpL5']
	)

l.bsc <- qnorm(p = 0.005, mean = Means.bsc, sd = SEs.bsc)
u.bsc <- qnorm(p = 0.995, mean = Means.bsc, sd = SEs.bsc)
ylim = c(-0.05, 0.31)

pdf('Figure6.pdf', width = 11, height = 4)
	plot(x = -2:5, y = Means.bsc, xlim = c(-2.5, 5.5), ylim = ylim, xaxt = "n", yaxt = "n", ylab = '', las = 2, xlab = '')
	lines(x = -2:5, y = Means.bsc, lty = 3)
	abline(h = 0, lty = 2)
	abline(v = 0, lty = 2)
	arrows(x0 = -2:5, y0 = l.bsc, y1 = u.bsc, angle = 90, code = 0, lwd=2)
	axis(side = 1, at = -2:5, labels = c('2 Years Prior', '1 Year Prior', 'Year of Event', '1 Year After', '2 Year After', '3 Year After', '4 Years After', '5 Years After'), cex.lab=4.2)
	axis(side = 2, at = seq(-0.1, 0.5, by = .1), labels = seq(-0.1, 0.5, by = .1), cex.lab=4.2, las = 3)
	title(ylab="Effect on State Repression", line=2, cex.lab=1.2)
	title(xlab="Time Passage Relative to Successful Coup", line=3, cex.lab=1.2)
dev.off()

